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Abstract 

We outline a general theory for the analysis of flow-distributed standing and travelling wave 
patterns in one-dimensional, open plug-flows of oscillatory chemical media. We treat both the 
amplitude and phase dynamics of small and large-amplitude waves, considering both travelling 
and stationary waves on an equal footing and emphasizing features that are generic to a variety 
of kinetic models. We begin with a linear stability analysis for constant and periodic boundary 
forcing, drawing attention to the implications for systems far from a Hopf bifurcation. Among 
other results, we show that for systems far from a Hopf bifurcation, the first absolutely unstable 
mode may be a stationary wave mode. We then introduce a non-linear formalism for studying 
both travelling and stationary waves and show that the wave forms and their amplitudes depend 
on a single reduced transport parameter. Our formalism sheds light on cases where neither the 
linearized analyis nor the kinematic theory of phase waves give an adequate description, and it can 
be applied to study some of the more complex types of bifurcations (Canards, period-doublings, 
etc.) in open flow systems. 

PACS numbers: 82.40. Ck, 82.40.Bj, 47.70.Fw 



INTRODUCTION 



Among the variety of known mechanisms of spontaneous pattern formation in spa- 
tially extended chemical systems, the flow distributed oscillation (FDO) mechanism 
* 16| differs from the Turing and related mechanisms 3] and from the differential flow 



instability H3 

in that it does not require anydifferential transport of the reacting species. 
The FDO mechanism, which was predicted physically interpreted and experimentally 
verined0-@, operates in open flow systems when the chemical kinetics is intrinsically oscil- 
latory and the inflow boundary condition plays an essential role. It is potentially relevant to 
biological pattern formation involving an axial growth, since an open flow system is related 
by a coordinate transformation to a linearly growing medium such as a plant stem or animal 
embryo. [Il|-ji3 

The governing equations of the one-dimensional open reactive flow studied are the 
reaction-diffusion-advection (RDA) equation 

<9U dU d 2 V 

— = f(lJ;C)-v— + D—. (1) 

together with the boundary condition U(0, t) at the inflow. Here XJ(x,t) is an N- 
dimensional vector of dynamical variables (concentrations of species), f (U;C) is the vector- 
valued rate function which depends on one or more control parameters C, v > is the flow 
velocity and D is the diffusion coefficient. In general D and v can be matrices, allowing 
for differential transport, but here we focus on the case without differential transport, so 
that v and D are real scalars. We take the length of the reactor to be L and impose no-flux 
boundary conditions at the outflow: dU / dx\ x=L = 0. We are interested in the case where 
f(U;C) has a stable limit cycle and at least one unstable fixed point. 

Insight into the physical mechanism of flow-distributed waves can be gained by consid- 
ering the kinematic limit [7] of zero diffusion. In this case equation (JTJ) can be written 
as 

f = f(U;C), (2, 
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where we have introduced the advective derivative d/dt = d/dt + vd/dx. Then the evolu- 
tion of an advected fluid element is the same as that of the well-mixed system: each fluid 
element is an independent oscillator or batch reactor, and the inflow boundary condition 
serves to establish the phase of these oscillators. [31] If the boundary condition is constant, 
for example, then all of the oscillators enter the reactor with the same initial phase and sta- 
tionary waves result from the recurrence of the same phase at equally spaced downstream 
positions. On the other hand, an oscillating boundary condition leads to upstream or 
downstream travelling waves, which have also been verified experimentally. [ll]-[16j While 
the kinematic limit is helpful for understanding the essential physics, there are significant 
deviations from its predictions when diffusion becomes significant D / 0. These deviations 
affect the wavelengths as well as the amplitudes and shapes of flow distributed waves. 
The path in phase space of a volume element moving through the reactor may differ from 
that of the well-stirred system with the same initial conditions. |H Sufficiently strong diffu- 
sion can extinguish the waves due to the diffusive mixing of adjacent fluid elements which 
are oscillating out of phase. 

An alternative aproach to this kinematic or phase dynamics picture is the original ^ 
linear stability analysis which views the flow-distributed oscillations as arising from growing 
perturbations of an unstable uniform stationary state U(x, t) = Uo where Uo is an unstable 
fixed point of the underlying kinetics governed by f (U;C). This linearized analysis is useful 
for predicting which small perturbations initially grow and which do not and it has been 
used in a variety of systems to determine thresholds for the formation of wave patterns 
(for example, the bifurcation from growing to evanescent stationary waves.) However, 
the linearized analysis is not sufficient to determine the behavior of the waves when their 
amplitude becomes large. As was pointed out in ref. the linearized analysis is most 
useful when the kinetic system is not far from a supercritical Hopf bifurcation. In this case, 
it can give a good approximation to the wavelength of the fully developed nonlinear waves, 
although by itself it is inadequate to predict their amplitude or shape. Chemical oscillators, 
however, have interesting dynamical regimes far from a Hopf bifurcation, showing markedly 
non-sinusoidal or relaxation oscillations. In such cases, the fixed point enclosed by the limit 
cycle may be an unstable node rather than a focus, in which case a linearized analysis near 
the fixed point does not reveal the intrinsic oscillatory behavior at all. The behavior of 
large perturbations can then be radically different from that of small ones. 

Our goal in this work is to develop a general formalism that describes both stationary and 
travelling flow-distributed waves including their nonlinear behavior and deviations from the 
kinematic limit. This formalism provides a unifying framework for understanding previous 
results and it can be applied to arbitrary kinetic models. In contrast to previous studies of 
particular kinetic models[H@i[I3, our emphasis is on generic phenomena. In particular 
we hope to shed light on systems with relaxation or other non-sinusoidal oscillations, where 
the unstable fixed point is a node rather than a focus. In Section II we discuss the linear 
stability analysis of a fixed point with an emphasis on generic features of the two types of 
fixed points, unstable foci and unstable nodes. We derive a general dispersion relation for 
small-amplitude disturbances and use it to analyze the response to time-dependent pertur- 
bations imposed at the inflow. We show that a band of frequencies is spatially amplified, 
and that this band becomes narrower and sharper with increasing diffusion. We derive 
general expressions for the thresholds separating absolute from convective instability as well 
as the threshold for extinction of stationary waves and show that the latter threshold dis- 
appears in the case of an unstable node. If the system is far from a Hopf bifurcation, we 
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show that, in contrast to previous examples, it is possible for sustained stationary waves to 
arise through an absolute instability from a transient perturbation. In Section III we use a 
travelling wave ansatz to introduce a reduced one- dimensional ordinary differential equation 
(ODE) that describes the fully nonlinear stationary and travelling waves. Numerical solu- 
tions of this equation are easier to obtain than those of the partial differential equation (JTJ) 
and they allow us to obtain the nonlinear dispersion relation for the fully developed, large 
amplitude waves. A key result is that the amplitude and shape of the wave depend only on 
the reduced transport parameter T = D/(v — c) 2 where c is the phase velocity of the wave. 
It is T that governs the strength of deviations from the kinematic limit. We show that in 
the case of relaxation oscillations these deviations are qualitatively different from those in 
the quasi-sinusoidal case, and we provide a physical interpretation of the dependence on 
T. We illustrate the application of our techniques by some numerical results. Finally, 
in Section IV we summarize our conclusions and describe the prospects for applying our 
techniques to more complex types of bifurcations (Canards, period doublings, etc.) in open 
flow systems. 

Throughout the paper we use as an illustrative model the van der Pol or FitzHugh- 
Nagumo (FN) oscillator Q EH @ 

f - < x ~ x3 ~ >•> < 3 > 

dY 

w-- Y + 10X ' 

which is described in more detail in Appendix A. The FN model exhibits the generic 
features we are interested in: it provides examples of quasi-sinusoidal oscillations changing 
to relaxation oscillations and an unstable focus changing to an unstable node as e is increased. 
Similar qualitative features occur in the Brusselator, Oregonator and other kinetic models. 
Appendix B describes the numerical approaches used and challenges encountered in solving 
the 1-D ODE described in Section III. 



II. GENERIC LINEARIZED ANALYSIS NEAR A FIXED POINT 

In this section we pursue the approach which views flow-distributed waves as arising 
from instabilities of the spatially uniform solution U(x,t) = Uo of equation ((H). The 
systems of interest possess an unstable fixed point and a stable limit cycle. The instability 



of the uniform state may be either convective or absolute. [23] 24] 25] We now consider small 
perturbations u(x, t) 

U(x,t) = U + u{x,t). (4) 
of the homogeneous fixed point Uo- In the linearized approximation, u(x,t) obeys 

- = mo )u-v- + D^ (5) 

where J is the Jacobian matrix <9f(U)/<9U. Let us denote the eigenvectors and eigenvalues 
of the Jacobian by ^ and Xj respectively (j G {1..N}) and expand the perturbation u(x, t) 
in the eigenbasis as 

N 

u(x,t) = ^2uj(x,t)^j. (6) 
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Substitution into the linearized equation gives a separate equation for each component: 

dt 33 dx + dx 2 ' [) 

We are interested in the dynamics along unstable eigenvectors, i.e., ones for which the 
real part of Xj is positive. Xj may be one of a pair of complex conjugate eigenvalues 
X± = a±i(3 with a positive, in which case U is an unstable focus, or it may be purely real 
and positive, as for an unstable node. In general, the former case occurs in the neighborhood 
of a supercritical Hopf bifurcation, while the latter may occur farther from the bifurcation. 
The component equation associated with the eigenvector ^ can be written as 

+ < 8 » 

Without loss of generality we can assume that (3^0. The case of a purely real eigenvalue 
is included by setting (3 = 0. We seek complex exponential solutions of this equation of the 
form 

u(x, t) = Aexp(iut + ikx). (9) 

(From here on, the subscript j is suppressed.) Both ui and k can be complex: uj = uj t + ioji 
and k = k r + iki. Our sign conventions are such that positive values of u>i or ki correspond 
to solutions that damp with time or downstream distance respectively. Substitution of the 
complex exponential ansatz @ into equation (|SJ) gives the dispersion relation 

iuo = a + i(3 - ivk - Dk 2 . (10) 

We now consider modes with Ui = 0, which represent the response to a sinusoidal forcing 
at the boundary with frequency uo r . With this restriction to real frequencies, the real and 
imaginary components of the dispersion relation (jl(Jj) read 

vk t + D(k 2 -k 2 r ) = -a, (11) 



vk r + 2Dk r ki = (3 - u T . (12) 

If u r — [3 (forcing at the natural oscillation frequency), then eq. (J12)) can be satisfied by 
setting k r = 0. Eq. (|TT|) is then reduced to a quadratic equation in k iy which always has 
a negative (spatially growing) solution provided D/v 2 < l/4a. (As we will discuss below, 
D/v 2 = l/4a marks the threshold between convective and absolute instability of the fixed 
point. Beyond this threshold the dispersion relation cannot be solved with a purely real 
to.) >>From this we learn that a perturbation with uj r = (3 always yields a growing mode 
(the Hopf mode) in which the whole reactor oscillates in phase. In fact, we shall see that 
(3 is the midpoint of a band of amplified frequencies. The edges of the amplified band can 
be derived by setting ki = in eqs. ()11|) and (|T^|) and solving for u r . The result is that 
sinusoidal perturbations are spatially amplified if 



av 2 „ / av 2 



P-)l^ r <"r<(3+ ] l^ r (13) 

and they are damped otherwise. If 

Da . . 

> W (14) 

then u) r = falls outside the range of amplified frequencies and the stationary waves decay. 
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A. Threshold between Absolute and Convective Instability 



An important distinction in the dynamics of open flow systems is the one between absolute 
and convective instability. If a state is convectively unstable, then perturbations of the state 
grow in a reference frame moving with the flow velocity, but for any fixed position in the 
laboratory frame the disturbance eventually decays. The effects of a perturbation of finite 
duration are eventually washed downstream and out of the system and the reactor returns 
to its initial unperturbed state. In this case, the initial conditions at t = do not influence 
the asymptotic late-time behavior, which is instead controlled by the upstream boundary 
condition. If a state is absolutely unstable, on the other hand, then localized perturbations 
grow with time at a fixed position in the laboratory frame, propagating upstream as well 
as downstream. This means that perturbations applied at t = continue to affect the state 
at late times. From an experimental point of view, the continuing effect of noise in the 
initial conditions makes an absolutely system more complicated. Since the disturbances 
propagate in both directions, the system resembles a pure reaction-diffusion system where 
the flow plays a subordinate role. The critical slowing-down near the convective/absolute 
instability threshold adds to the experimental difficulty. A convectively unstable system 
can be controlled more readily by manipulating the inflow boundary condition. 

The distinction between the two cases can be expressed in terms of the group velocity 
of propagating disturbances. (H [HI [HI If there are modes with positive growth rate and 
upstream (negative) group velocity then the instability is absolute — otherwise, all growing 
disturbances are washed downstream. The threshold between the two cases occurs when a 
mode with zero group velocity dui/dk = has exactly zero growth with respect to time, 
Ui = 0. This represents a persistent disturbance that is not washed downstream. From 
the dispersion relation (jl(Jj) the condition of zero group velocity is 

duj 

= — = -w - 2Dk = 
dk 

or 

k = w- (15) 

giving a purely imaginary wave number which can be substituted into both components of 
the dispersion relation (fTTj) and ()12j) . The threshold for absolute instability occurs when 
the temporal growth rate — Ui is exactly zero for this zero-group- velocity mode, so from the 
real component (fTTj) we get 

2 2 
V . V , 2 V 



= = a h D( — Y = a 

2D K 2D J AD 



or 



(16) 



D 1 

v 2 Aa 

From the imaginary component (JT2jl we get that u r = f3, which tells us that the first mode 
to go absolutely unstable is the Hopf mode. 

Note that eq. (|16|) is the same threshold mentioned above in discussing the solutions of 
(jll|) and (|12j) . It is the threshold beyond which the most unstable mode grows with time 
as well as space and therefore its frequency cannot be purely real. 
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FIG. 1: General solutions to the dispersion relation for small-amplitude oscillations near an unsta- 
ble fixed point. All frequencies are scaled relative to the inverse growth time scale a. Recall that 
negative values of the imaginary wavenumber Qi represent growing modes. 

B. Universal solutions of the dispersion relation 

To illustrate the above results generically, it is helpful to express the dispersion relation 
in terms of dimensionless variables. This requires a choice of a characteristic time scale. It 
is convenient to use the growth time scale t = 1/a and introduce the variables Q = u/a, 
Q = Q r + iQi = vk/a, e = Da/v 2 and Q = [3 /a. For real Q the two components (fTTj) and 
(fT2Jl of the dispersion relation can be rewritten as 

-Qi + e(Q 2 r -Q 2 ) = l, (17) 

Q r + 2eQ r Qi = Q-Q . (18) 

Numerical solutions for Q r and Qi as functions of Q — Qq ar e plotted in figure for several 
values of e. The family of curves illustrates the qualitative features of the response to 
sinusoidal forcing: the band of amplified frequencies becomes becomes sharper and narrower 
with increasing e, and the plot of the wavenumber Q r versus Q becomes more strongly 
nonlinear. The maximum growth rate always occurs at Q = Qq. The onset of absolute 
instability occurs when e — 1/4 and is marked by the appearance of a cusp in Qi and a 
vertical jump in Q r at Q = Qq. In the kinematic limit e — > all frequencies are amplified 
equally and Q r is a linear function of Q. 

C. Comments on the unstable node case 

In the preceding paragraphs we have derived results for the linear stability analysis of 
perturbations along a single eigendirection near a fixed point of the dynamics. The results 
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are formally independent of whether the fixed point is a focus or a node. The latter case can 
be described simply by considering (3 = 0. However, there are important qualitative differ- 
ences between the two cases. In the case of an unstable focus, the eigenvalues are complex 
conjugates and trajectories spiral away from the fixed point. There is a single growth rate 
a. For a node, on the other hand, there is no intrinsic oscillatory behavior near the fixed 
point and any oscillation must be imposed externally or arise from nonlinear effects at larger 
amplitudes. There are generally two eigendirections associated with distinct growth rates 
ctj and therefore different dispersion relations and instability thresholds for perturbations 
along the two eigendirections. Perturbations grow most rapidly along the most unstable 
eigendirection and therefore this eigendirection is physically the most important. 

Another difference between the two cases is that for a focus ((3 7^ 0) the most unstable 
mode is always the Hopf mode which is a uniform oscillation at frequency (3. The band 
of amplified perturbations is centered on this nonzero frequency, which is also the first 
mode to become absolutely unstable. When (3 = 0, on the other hand, the most unstable 
mode (and the first to become absolutely unstable) is a stationary mode with u r — 0. The 
linearized dispersion relation predicts a real component of the wavenumber k r = for this 
mode. However, when the amplitude grows, the nonlinear terms grow rapidly with the 
result that large-amplitude oscillation and stationary waves of finite wavelength take over. 
The different nature of the absolute instability threshold between the two cases is illustrated 
with numerical results in figure El 1 The emergence of sustained stationary waves due to 
absolute instability from a time-limited perturbation (fig. Eb) is i n sharp contrast to previous 
examples[§]0 in which the first absolutely unstable mode was time-dependent (as in fig. 
12k) and stationary waves in the absolutely unstable regime only occured in a thin boundary 
region, if at all, and only in response to a steady boundary condition held away from the 
fixed point. 



D. Bifurcation Loci in Control Parameter Space 

In a reactive flow experiment, it is typically the flow velocity which is under the exper- 
imenter's control while the diffusion constant is fixed. One focus of previous studies has 
therefore been the determination of critical flow velocities at which the FDO patterns change 
qualitatively. Two important thresholds are the threshold for the formation of undamped 
stationary patterns vt and the threshold between absolute and convective instability, v c , 
using the notation of ref. ■ General expressions for these threshold velocities are readily 
obtained from (j!4j) and (fTTtj) : 

MV-JWi (19) 

y/a{C) VD 

v c (C) = 2^(Cj-±=, (20) 



where the dependence on kinetic control parameters C has been made explicit. These 
thresholds were plotted as functions of a control parameter for particular systems in refs. 



1 Numerical solutions of the RDA equation were obtained using a simple first-order discretization. The 
time and space grids were adjusted according to the characteristic time and space scales of the system 
being studied. For a sufficiently fine grid, the results were verified to be insensitive to the grid size. 
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FIG. 2: Onset of absolute instability in the FN system. The gray levels are proportional to 
X(x,t). The inflow boundary conditions are X(0,t) = —0.001, < t < T, X(0,t) = otherwise, 
where T = 2 in case (A), 0.1 in case (B). If the system were convectively unstable, the effects of 
the perturbation would be washed out of the system, but instead both systems are just above the 
threshold of absolute instability and the effects of the disturbance persist. The pulse perturbation 
has a fourier spectrum of frequencies, of which the most unstable grows most rapidly. In case 
(A) (e = 2, v = 0.1, D = 0.0052) the fixed point is a focus and the first absolutely unstable 
mode is oscillatory. In case (B) (e = 50, v = 1, D = 0.01) the fixed point is a node and the first 
absolutely unstable mode is a stationary wave. Although the linearized theory correctly predicts 
that the persistent mode is a stationary wave, it does not predict the finite wavelength, nor does 
it adequately describe the downstream front of the disturbance. 



0,0] and We can now understand some generic features of the shapes of these bi- 
furcation loci observed in those works. Near a supercritical Hopf bifurcation, vt diverges 
and v c —>■ because a — » 0. In a range above the Hopf bifurcation, then, vt > v c and 
there are three dynamical regimes as the flow velocity is varied. At high velocities, there 
are sustained stationary waves which can be excited by fixed boundary conditions. For 
v c < v < vt, the stationary waves are evanescent and penetrate only a limited distance into 
the medium. After any initial transients 2 are convectively washed out of the system, the 
remainder of the system returns to the fixed point. Thus there is a range of v values for 
which the unstable fixed point is effectively stabilized in the presence of constant boundary 
conditions. (However, oscillatory perturbations in the amplified band still result in un- 



2 These transients can themselves have complex and interesting structures which we do not consider here; 
see, e.g., ref. 10]. 
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damped travelling waves.) As a control parameter is tuned farther away from the Hopf 
bifurcation, a(C) generally increases, and the two curves vt{C) and v c (C) may cross (see, 
for example, figure 1 of ref. [2J) so that for some values of C, vt < v c . In this case there is no 
intermediate regime of evanescent waves followed by return to the fixed point as described 
above. Instead, as the velocity is lowered, the absolute instability threshold is reached first. 
Below v c , the effects of initial conditions at t = are not washed out of the system and the 
late-time behavior may not be controllable by the boundary condition. The most unstable 
modes dominate, and stationary waves may occur only in a boundary layer with constant 
boundary conditions. It is possible that the crossing of the two thresholds (vt < Vc) ac- 
counts in part for the failure to observe evanescent stationary waves in experiments. 0]]^ In 
some cases, as C is tuned still further from a Hopf bifurcation, /3(C) may decrease to zero, 
and the threshold vt vanishes as the fixed point becomes an unstable node. In this case, 
stationary waves may be triggered by an absolute instability as mentioned above. 



III. NONLINEAR, LARGE AMPLITUDE WAVES 

In this section we derive a reduced ordinary differential equation (ODE) that describes 
both stationary and travelling wave solutions of the RDA equation (|1J) and that applies to 
situations where neither the linear stability analysis nor the kinematic limit give an adequate 
description. We show that the amplitude and wave form depend on a single reduced 
transport parameter which characterizes the degree of departure from the kinematic limit. 
The application of our formalism is then illustrated by some numerical examples. As in 
the linearized analysis, we find that the two cases near and far from a Hopf bifurcation give 
qualitatively different wave behavior. 

We make the ansatz u(x,t) = u c (x — ct) which is analogous to the D'Alembert solution 
of the wave equation. It represents a generic wave (not necessarily periodic) travelling 
downstream with velocity c and depending only on the combination ( = x — ct. A growing 
or decaying travelling wave is not strictly described by this ansatz unless c = (since a 
spatially changing amplitude implies a dependence on x and not purely on Q but we expect 
that it describes the asymptotic behavior of such a wave when the amplitude saturates. 
Substituting this into the reaction-diffusion-advection equation (0) gives a one-dimensional 
ODE: 

= f(u)-{v-c)^ + D^. (21) 
With a change of variable = v — c) we obtain 

„, . du ^d 2 u 

° = f(u >-5c + r 3P (22) 

where T = D/(v — c) 2 . T represents the effective strength of the diffusion term for a given 
wave. If (J22|) has a periodic solution u(£') = u(£' + A), then the period A in terms of (' is 
related to the frequency and wave number of the corresponding travelling wave by 

t c 27T 

u = ck = — — , 23 

c — v A(l ) 

where we have made explicit the dependence of A on V. 
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Some limiting cases are instructive and help to furnish a physical interpretation of the 
particular combination of parameters T. First, note that the kinematic limit r — > can 
be approached in several ways, either by decreasing D or by increasing \v — c\ . The limits 
c — > ±00 (infinite phase speed) both correspond to oscillations which are uniform in space. 
Diffusion can have no effect on a system with no spatial gradients, and so it is understandable 
that the effective diffusion vanishes in this limit. On the other hand, as c approaches v 
from either direction T diverges. This divergence can be understood qualitatively if one 
considers the diffusionless limit as a "zeroth order" approximation to the phase dynamics. 
In this limit, A = A(0) is a constant and the physical wavelength 27r//c shrinks to zero as 
c — > v. If a nonzero diffusion is then turned on, then it is natural to expect that it creates 
more mixing between crests and troughs when the waves are close together than when they 
are far apart. As the wavelength shrinks, the effect of diffusion on the waveform increases. 

The special case c = describes stationary waves, which result from a steady-state 
boundary condition. The combination of parameters controlling the behavior of stationary 
waves is the same combination identified in the linearized analyisis, namely r| c=0 = D/v 2 . 
In the case of complex eigenvalues, the condition for the existence of undamped periodic 
stationary waves is the one already derived using the linearized analysis: T < a/ (3 2 . Above 
this threshold all solutions spiral into the fixed point. However, since the nature of solutions 
to (|22p is governed only by T, one can obtain a more general threshold for undamped periodic 
travelling waves with phase velocity c. Such waves are undamped only if 



For given values of D and of the kinetic control parameters that determine a and /3, there is 
an interval of phase velocities, surrounding the flow velocity, within which travelling waves do 
not propagate. For the case of an unstable node or (3 = 0, on the other hand, the threshold 
diverges and so it appears that there is no such excluded band of phase velocities. We 
will further examine this issue below. 

A. Numerical solutions and qualitative features of nonlinear waveforms 

The equation ()22j) is useful for several reasons. First, it shows that the essential 
features of travelling and stationary waves (the amplitude and the waveform) depend on a 
single combination of transport parameters, F , and therefore a family of different waves are 
described by a single universal function. Second, as a general rule, ODE's can be solved 
with less computational effort than PDE's, and (}2*2"j) allows one to derive wave solutions 
without solving the full PDE (JTJ. Solution methods are described in Appendix B. (|2*2*jl is 
solved on a finite interval with boundary conditions, and the solution will in general include 
transients near the boundaries. In the case c = 0, these may represent actual transients 
of the stationary wave solution near the physical boundaries of the reactor. On the other 
hand, if we wish to describe the asymptotic behavior of travelling wave solutions with 
the transients should be ignored. 

Some examples of numerical solutions for the FN system are shown in figures El and 
I3J As expected, increasing T has the general effect of increasing the deviations from the 
kinematic limit. However, the behaviour differs qualitatively between the quasi-sinusoidal 
and the relaxation oscillation cases. In the former case, the phase space orbit remains 
approximately elliptical. Its period remains approximately constant while its amplitude 
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(24) 
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FIG. 3: Numerical solutions of equation (|22jl for the FN system with e = 2. The left column 
shows Y(() vs. -X'(C)) the right shows X(£) vs. £. The dotted and lighter-colored curves show 
solutions in the kinematic limit (or T = 0). The boundary conditions are fixed at a point on the 
local limit cycle at one end, free at the other end. For T < 0.28 ((A) and (B)) there is a periodic 
quasi-sinusoidal waveform. For T > 0.28 (C) the solution spirals into the origin. The period 
changes little as T increases. 

shrinks uniformly until, at the critical threshold T = a//3 2 , it vanishes into the fixed point. 
In the relaxation oscillation case, on the other hand, the limit cycle does not shrink to a 
point. Instead, as T increases, the period A increases apparently without bound, and the 
limit cycle remains elongated in the more unstable eigendirection while narrowing in the 
transverse direction. 

B. "Non-linear dispersion relation" 

By finding A(r) numerically for a range of values of values of T and then using the relation 
(I23|) . one can find the non-linear dispersion relation between the frequency u (which can 
be set by the forcing frequency of a perturbation at the inflow) and the phase velocity c 
of large-amplitude travelling waves at given values of the transport parameters D,v. In 
the quasi-sinusoidal dynamical regime near the Hopf bifurcation, A is nearly constant and 
approximately equal to the small-amplitude oscillation period 2tt//3. The uj-c relation 
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FIG. 4: Solutions of equation (|22|) as in figure 01 but for e = 50. As T increases, the period 
lengthens compared to the kinematic limit and the waveforms trace narrower loops in phase space, 
but they do not spiral into the origin at any finite T. 



is then approximately the same as that predicted in the kinematic limit, except that it is 
truncated at the cutoff frequencies (3 ± a/ 'av 2 / D where the amplitude vanishes (see figure 
|U|3) . Outside of this interval of frequencies (which is the same as the range of amplified 
frequencies in the linearized theory) there are only evanescent waves. Because of the inverse 
relation between u and v — c, there is, as mentioned above, a range of excluded phase 
velocities near v for which undamped travelling waves do not exist. 

In the relaxation regime, on the other hand, A varies quite strongly with T. In fact, 
numerical results suggest that for asymptotically large T it increases approximately linearly 
(see figure EI). This means that the frequency does not become infinite as c — > v but that 
instead there is a maximum. Such a maximum is seen in figure EK, which shows frequency 
versus phase speed for the case e = 50, D = 0.003, v = 1, based on the numerical data 
for A(r). This relation is radically different from the one for small-amplitude, linearized 
waves. The maximum frequency appears to be lower than the cutoff frequency obtained 
from the linear stability analysis. There is thus a range of frequencies for which the linear 
theory predicts a growing mode, yet there are no large-amplitude solutions described by 
corresponding to these frequncies. Numerical results described below suggest that within 
this frequency gap the small-amplitude waves are subject to a secondary instability and do 
not penetrate far beyond the inflow boundary. 
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FIG. 5: Scaled wavelength A as a function of T obtained from numerical solutions of the one- 
dimensional equation. A appears to increase approximately linearly for large T. The numerical 
results become more uncertain at longer wavelengths due to the finite interval of the solution (see 
Appendix B). 
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FIG. 6: (A) Phase velocity vs. frequency uj for the case e = 50, v = 1, D = 0.003 based on the 
numerical data of figure |SJ There is evidently a maximum frequency. Perturbations near this 
frequency generate waves with a speed close to the flow velocity. The phase velocity for linearized 
small- amplitude waves along the most unstable eigendirection is shown for comparison. Note that 
the cutoff frequency for small-amplitude waves occurs precisely when their phase velocity reaches 
1, as can be deduced from the dispersion relation. The maximum frequency for large-amplitude 
waves (approximately 90) appears to be lower than the linear cutoff frequency of 110. In the gap 
between these two frequencies, the linearized analysis predicts a growing mode but eq. (|22jl gives 
no solution with the correct frequency. (B) A similar plot for e = 2, v = 1, D = 0.0025. In this 
case the phase velocity for nonlinear waves is very close to the linearized prediction, and becomes 
closer as the cutoff frequency of 18.36 is approached. 
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FIG. 7: Growing travelling waves for e = 2, v = 1, D = 0.0025 generated by a sinusoidal boundary 
perturbation with u = 15. (A) X(x,t) represented as a gray level. (B) X(x) at t = 20. (C) 
Y{x) vs. X(x) at t = 20. The limit cycle of the underlying oscillator is shown as a lighter curve, 
and the cubic nullcline (thin dashed line) is also shown for reference. The actual waveform has 
a slightly smaller amplitude than the limit cycle of the well-stirred system. The saturated waves 
have a phase velocity of 1.65, giving T ~ 0.006. 



C. Evolution and asymptotic waveforms of growing modes 

In this section we present a few numerical solutions of the PDE (0) with sinusoidal bound- 
ary forcing (X(0,t),Y(0,t)) = (a cos cut, a sin tut) where a is a small perturbation amplitude 
(a = 0.05 for most examples). These examples illustrate the principles discussed above, 
including the role of the reduced parameter T and the growth of small perturbations into 
nonlinear waves that obey ()22|) . First, we illustrate the "equivalence" of waves with the 
same value of T. Figures [7| and |S] are space-time diagrams of the FN system at e = 2. The 
diffusion constants are different, and in one case to — 15 while to = in the other. But in 
both cases T « 0.006 and so the functions XJ(x) for fixed t have approximately the same 
shape and saturate at almost the same amplitude (slightly below the kinematic limit). 

Next, we examine numerical results for e = 50, v = 1 and D = 0.003, where the unstable 
fixed point is a node with two positive real eigenvalues. The predicted uo-c relation for these 
parameter values is shown in figure The boundary perturbation has components of 
equal size along both eigenvectors, but the component along the more unstable eigenvector 
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FIG. 8: Stationary waves in the FN system with e = 10, v = 1 and D = 0.006 have the same value 
of — as the travelling waves in figure [7| Their spatial waveform is similar in shape and amplitude 
and deviates from the kinematic limit (lighter curve) by the same amount. 

of course grows faster. Figure IH1 and ITUl show waves generated by perturbations with uj = 50 
and 80, respectively. In both cases, the perturbation grows to a nonlinear travelling wave 
that deviates significantly from the kinematic limit (with more pronounced deviations in the 
uj = 80 case.) The waveforms resemble those of figure EJ Figure ITT1 shows the results of a 
perturbation with uj = 85 which is near the maximum frequency for nonlinear waves in figure 
El In this case, the perturbation initially grows along the more unstable eigendirection, 
but it penetrates only a small distance into the medium before breaking up. Similar results 
were found for perturbations between uj ~ 85 and the linear cutoff frequency of uj ~ 110 
(For uj > 110 the perturbations are immediately damped.) Evidently the waves within this 
frequency range are subject to a secondary instability. The patterns which arise after the 
high-frequency waves break up appear similar to the pulsating waves observed in refs. 0] 
and The behavior of perturbations near the cutoff frequency in the case of an unstable 
node warrants further study. 
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FIG. 9: Waves generated by a sinusoidal perturbation with u = 50 in the FN system with e = 50, 
v = 1 and D = 0.003. The initial transient response to the switching on of the perturbation at 
t = is followed by a steady travelling wave with phase velocity c ~ 1.6, a value consistent with 
the oj-c relation plotted in figure The corresponding value of V is 0.008 (C) shows deviations 
from the kinematic limit — the shape of the local orbit is shown as the broken line. Compare (b) 
and (c) to figure E^. 



IV. CONCLUSIONS AND DISCUSSION 

We have attempted to give a general framework for understanding the behavior of flow- 
distributed waves in one-dimensional open flows of oscillatory media without differential 
transport, aiming at generic results. From the point of view of linear stability analysis, we 
examined some changes that occur generically when an unstable fixed point changes from 
a focus to a node, as typically occurs at sufficient distance from a Hopf bifurcation. We 
showed that when the fixed point is a node, the first mode to become absolutely unstable 
is a stationary wave mode. In this case, unlike examples previously noted, a stationary 
wave can be established by a perturbation of finite duration in time, at least if the system 
is not far above the absolute instability threshold. We derived general expressions ( eq. 
(113)0 for the band of frequencies that result in growing waves. The center of this band 
occurs at zero frequency if the unstable fixed point is a node, and at a non-zero frequency 
if it is a focus. The expression for the growing frequency band allows easy calculation of 
the threshold for the extinction of stationary waves. We then showed that the nonlinear 
behavior of periodic travelling or stationary waves reduces to a one-dimensional ODE (j22j) 
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FIG. 10: Waves generated by a perturbation with u = 80 show a larger deviation from the kinematic 
limit. The waves have c = 1.27 and T = 0.04. Compare (b) and (c) with figure HJx 



governed by the single parameter T — D/(v — c) 2 , which has a physical interpretation as the 
strength of diffusive mixing between peaks and troughs of a particular wave. The ODE can 
be solved numerically to derive the waveforms and obtain a relation between the frequency of 
driving at the boundary and the wavenumber and/or phase velocity of the waves generated 
by the perturbation. We examined deviations of waveforms from the kinematic limit, noting 
qualitative differences between the quasi-harmonic and relaxation cases. 

We illustrated our formalism by applying it to the FitzHugh-Nagumo toy model, but the 
tools we developed here can be applied to other kinetic models. In particular, we plan 
to use the current results to analyze FDO patterns of complex and chaotic oscillators in a 
future publication. The reduced one-dimensional equation can be applied to systems with 
multiple fixed points, subcritical Hopf bifurcations, Canards and bistable behavior, or other 
situations in which the linearized analysis gives only limited insight. 

Some other questions have been left open. The behavior of travelling waves near the 
cutoff frequencies in the case of an unstable node may be a fruitful subject for further study. 
More generally, we have only hinted at the possibility of secondary instabilities that may 
affect FDO waves. Also, in this work we have only considered regular travelling waves, not 
the pulsating waves observed in refs. [H and , although some of our numerical results (see 
figure ITT|) appeared to show pulsating waves arising from a secondary instability. 
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FIG. 11: Boundary perturbation frequency u = 85. A high-frequency travelling wave, oscillating 
almost entirely along the most unstable eigendirection, penetrates a limited distance into the 
medium before giving way to a pattern of much longer, irregularly moving waves. 



V. APPENDIX A: THE FITZHUGH-NAGUMO MODEL, QUASI-SINUSOIDAL 
AND RELAXATION OSCILLATIONS 

The FitzHugh-Nagumo (FN) model, also known as the van der Pol oscillator, is denned 

""3. 



by the following pair of equations 18] [21] [2 

j y 

-^iX-Xi-Y) (25, 
dY 

= -Y + aX + b. 

at 

It is not meant to be a realistic model of any chemical system, since its state variables include 
negative values, but it serves as a useful toy model that exhibits many generic features seen 
in real chemical systems, including bistability, excitability and ocillations of quasi-sinusoidal 
as well as relaxational character. 

The nullclines (Fig. n"2*|) are a cubic and a straight line, e is the ratio of time scales 
of motion along X and Y, a is the slope of the Y nullcline and b is its intercept with the 



3 This version of the model was used in ref. [l^ . 
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F-axis. relaxation oscillations occur when e is large. The number and location of the fixed 
points depends on the Y nullcline, i.e. on values of a and b. There may be either one 
fixed point or three, as shown in fig. A fixed point (X*,Y#) is always stable when it 

lies outside the extrema of the cubic, i.e. for X* < —l/y/3 and X* > l/y/3. Depending on 
the parameters (e, a, b) a single fixed point may be unstable and give rise to a limit cycle if 
-l/y/3 < X* < l/y/3. O ne case that has often been studied is that of a stiff system with 
a large time scale separation e ^> 1 and a single fixed point. In this case, as b is adjusted 
so as to move the fixed point closer to the origin, there is first a Hopf bifurcation at a fixed 
point near one extremum of the cubic nullcline, and then a Canard transition [26| in which 
the limit cycle changes rapidly from a small-amplitude one to a large-amplitude relaxation 
oscillation. However, we do not consider the Canard transition in this paper but instead 
follow reference in setting 6 = 0, thus ensuring a fixed point at the origin. By setting 
a = 10 we also ensure that there is only one fixed point, thus excluding bistable or excitable 
behavior. In the examples we discuss, we vary only e. 



(A) (B) 




— 1 O 1 " — 1 O 1 

x x 

FIG. 12: Nullclines for the FHN model for b = and (A) a > 1; (B) a < 1. In the latter case, 
there are three fixed points. 



The Jacobian evaluated at the origin is given by 

(26) 



e — e 
a — 1 



and its two eigenvalues are given by 



A ± = ^±^v / (l + e) 2 -4ea. (27) 

The eigenvalues are complex if and only if 

(1 + e) 2 1 1 e 
a> —^ = 4-e + 2 + 4- 

Otherwise, both eigenvalues are real, and their signs are the same if a > 1. If the eigenvalues 
are complex, then their real part is positive when e > 1. When a = 10, the eigenvalues 
are real and positive for all e > e CT n ~ 38. Figure UBI shows the real and imaginary parts 
of the two eigenvalues A± together with the angular frequency ujlc = 2n/T of the stable 



20 



limit cycle which exists for all e > 1. The Hopf bifurcation occurs at e = 1 where Re(A±) 
becomes positive. In the immediate vicinity of the Hopf bifurcation, the frequency of the 
limit cycle is identical to the imaginary part Im(A+). At e cr u, however, Im(A + ) vanishes and 
the eigenvalues become real. They are degenerate at the critical point but quickly become 
different as e increases further. Roughly speaking, the two real eigenvalues above e crit 
correspond to two different inverse time scales: a slower one for motion in the Y direction 
and a faster one for motion in the X direction. It is this separation of time scales that 
distinguishes relaxation oscillations from quasi-sinusoidal ones. The frequency of relaxation 
oscillations is determined primarily by the slower of the two time scales. As e increases, the 
qualitative character of the oscillations changes from approximately sinusoidal to relaxation 
oscillations, as shown in figure El Figure ITU shows trajectories starting from points close to 
the fixed point, therefore it also illustrates visually the change in character of the fixed point 
from an unstable focus (trajectories spiral away from the fixed point) to a node (trajectories 
leave along the most unstable eigenvector). Note that although there is a sudden change in 
the eigenvalues and eigenvectors near the fixed point at e cr i t , the associated change in the 
nonlinear limit cycle is gradual. 




FIG. 13: Inverse time scales as functions of e for a = 10. Dotted line: Im(A+); solid lines: 
Re(A±); dashed line: 2tt/T for the limit cycle. At the critical value e c ~ 38, the eigenvalues 
become real. 



VI. APPENDIX B: NUMERICAL SOLUTIONS OF THE REDUCED ONE- 
DIMENSIONAL EQUATION (USD 

Here we discuss the solution of the reduced ODE (|22|) : 

o = f(u) - — + r c — 



In the kinematic limit r — > 0, this equation reduces to a first-order equation, identical in 
form to that of the dynamics of the well-stirred system. Solutions of this first-order equation 
with general initial conditions approach the stable limit cycle of the well-stirred system. 
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X 

FIG. 14: Trajectories of an initial condition close to the origin, for a = 10 and three different 
values of e. (A) e = 1.3 is close to the Hopf bifurcation and the oscillations are quasi-sinusoidal. 
(B) e = 20 is an intermediate value, between the Hopf bifurcation and e c ~ 38. (C) e = 100 is 
well above e c and the limit cycle has the character of a relaxation oscillation. The change in the 
shape of the limit cycle is a gradual one, in spite of the rapid transition in the eigenvalues at the 
fixed point. 

However, for any nonzero value of T the equation is second-order and anti-dissipative 4 , 
so that the initial-value problem leads to unbounded solutions for most choices of initial 
conditions u(0) and u'(0). To exclude these unphysical solutions, a boundary condition 
must be imposed. Thus, the equation should be solved as a boundary-value problem on a 
finite interval < ( < L. Our procedure was to impose a fixed boundary condition on the 
left, u(0) = Ui n and a free boundary condition on the right, u'(0) = 0. In the case c = 
the boundaries correspond directly to the the physical boundaries of the plug-flow reactor. 

For the examples studied, we found that for moderate values of T and for sufficiently large 
L the solutions of the boundary problem behave qualitatively like solutions of a first-order 
initial value problem with an attractor. In other words, after some transient behavior at 
small ( which depends strongly on u in , the solutions settle either to a fixed value or to 
a periodic behavior with an intrinsic period A, 5 which depends on F but does not depend 



4 If the equation is rearranged to isolate the second-order on one side, the analogy with the equation of 
motion of a point particle shows that the "force" contains an anti-damping term, and the term — f(u) is 
also of the "wrong" sign, tending to push the particle away from the stable limit cycle of the well-mixed 
system. 

5 A chaotic attractor is also possible. We plan to discuss this in a future publication. 
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sensitively on Uj n or on L. The free boundary condition on the right affects the solution 
only in a small interval near the right boundary, i.e., the boundary can be moved to a 
larger value of ( without changing the solution on most of the interval < ( < L. In the 
c = case, the entire solution, including the boundary transients, is physically meaningful 
as part of a stationary wave pattern in the reactor. For c 7^ 0, boundary conditions at fixed 
( values are not directly equivalent to boundary conditions at fixed x, but the attractors 
reached by the solutions can be interpreted as the asymptotic shapes of travelling waves in 
the medium. 

In order to obtain numerical solutions of the boundary value problem, we used a col- 



location algorithm included in the Matlab software package. [27| This algorithm requires 
an initial guess for the solution which is then adjusted to satisfy the differential equations 
and the boundary conditions within a specified tolerance. For long solution intervals (large 
multiples of A) the algorithm may fail to converge unless the initial guess is close to the 
actual solution. We used two procedures for iteratively obtaining a solution: 

1. Use a solution of the initial- value problem of the first-order, T = system as an initial 
guess for a relatively small value of T. Then increase V iteratively to the desired value, 
using each solution as the guess for the next value of T. This procedure was used, for 
example, to generate the solutions at a range of values of T in figure El 

2. Sometimes it is more convenient to approximate the solution with piecewise solutions 
on a series of overlapping intervals. The procedure is as follows: First, solve the 
boundary value problem with the desired boundary condition Uj n at ( = and free 
boundary conditions at a relatively small L which is neither too much larger nor too 
much shorter than A. Obtain a solution Ui(£) on that interval. Then evaluate that 
solution at ( = L/2 and use ui(L/2) as the boundary condition for a new solution 
on the interval L/2 < ( < 3L/2. Continue this procedure on a series of overlapping 
intervals. If L is not too large, then the initial trial solutions need not be close to 
the final ones, and if L is not too small they are not sensitive to the free boundary 
condition at the right, so that the overlapping solutions should be approximately the 
same except very near the boundaries. Stitched together, the piecewise solutions 
approximate a solution on a longer interval. 

Large values of T (where large means significantly larger than l/4a, the threshold of 
absolute instability in the c = case) often presented computational challenges because the 
solutions become more sensitive to the right boundary condition, and large intervals were 
needed in order for an attractor to appear. This is the source of some of the numerical 
jitter in the data in figure El 
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